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We show that an attempt to compute numerically a viscous flow in a domain with a piece-wise smooth 
boundary by straightforwardly applying well-tested numerical algorithms (and numerical codes based on their 
use, such as COMSOL Multiphysics) can lead to spurious multivaluedness and nonintegrable singularities in the 
distribution of the fluid's pressure. The origin of this difficulty is that, near a corner formed by smooth parts 
of the piece-wise smooth boundary, in addition to the solution of the inhomogcncous problem, there is also an 
eigensolution. For obtuse corner angles this cigcnsolution (a) becomes dominant and (b) has a singular radial 
derivative of velocity at the corner. A method is developed that uses the knowledge about the eigensolution to 
remove multivaluedness and nonintegrability of the pressure. The method is first explained in the simple case of 
a Stokes flow in a corner region and then generalised for the full-scale unsteady Navier-Stokes flow in a domain 
with a free surface. 



1. Introduction 

In many problems of computational mechanics it is necessary to consider domains with piecewise smooth 
boundaries, in particular the situation where two parts of the boundary locally form a two-dimensional wedge. 
In elasticity such problems occur in fracture mechanics when one considers how a crack propagates into a 
material (see Moes et al. |1999 1 . A typical example from fluid dynamics is found in the process of dynamic 
wetting (see Dussan|1979 Blake 2006 1, in which a liquid spreads over the surface of a solid, with the liquid- fluid 
free surface forming what is referred to as a 'contact angle' with the solid substrate. The mathematical problems 
in wedge regions that one encounters in fluid and solid mechanics have much in common; below we shall consider 
some generic computational issues in the context of fluid mechanics. 

It has been known for a long time (see |Huh fc Scriven|197l" I that classical boundary conditions, when applied 
in dynamic wetting, i.e. in the situation where the 'contact line' formed by the free surface with the solid 
boundary moves with respect to the solid, lead to a physically unacceptable outcome, namely a nonintegrably 
infinite force acting between the liquid and the solid. Termed the 'moving contact- line problem', this difficulty 



has become the subject of many theoretical works (see Ch. 3 of Shikhmurzaev (20071 for a review). 

The common feature of almost all theoretical models proposed for the moving contact-line problem is that 
the no-slip boundary condition on the solid surface, i.e. the condition for the tangential component of the fluid's 
velocity, is replaced by a more general condition formulated for the tangential stress (and hence allowing for 
slip) . The impermeability condition for the normal component of velocity on the solid surface as well as both the 
zero tangential- stress and impermeability conditions on the free surface remain intact. As a result, one arrives 
at a mathematical problem, locally in a wedge region, with different boundary conditions on the two sides, 
both of which involve first-order differential operators applied to the bulk velocity component tangential to the 
boundary and linear homogeneous algebraic conditions for the normal component. 

In the present work, we show that an attempt to directly apply well-tested standard numerical methods to 
this problem can lead to a fundamentally unsatisfactory outcome. As is shown, even in the simplest case of a 
steady Stokes flow in a corner region, irrespective of a particular numerical implementation, in a certain range 
of wedge angles the computed bulk pressure is both nonintegrably singular and multivalued as one approaches 
the corner. This outcome has a serious 'practical' implication: it becomes impossible to compute a numerical 
solution to the corresponding free-surface problems with moving contact lines, as no capillary pressure would 
be able to balance a nonintegrable bulk pressure. Alternatively, if the pressure singularity and multivaluedness 
are suppressed using numerical means, the resulting 'solution' becomes mesh-dependent and hence is by no 
means a uniformly valid approximation of the solution of the original partial differential equations. Neither of 
these options is satisfactory, and the persistent nature of the problem makes it necessary to handle it robustly, 
beginning by finding out its origin. 

The structure of the present paper is as follows. In Section [2] we formulate the problem to be solved in the 
simplest case of a Stokes flow in a corner region and in Section [3]describe the local analytic asymptotic behaviour 
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of the solution near the corner. In Section [4] numerical results are presented. They show that the numerical 
solution is in excellent agreement with the analytic asymptotics for acute corner angles and completely incorrect 
for obtuse angles. In particular, the fluid's pressure appears to be multivalued at the corner with a noninte- 
grably singular behaviour near it. As demonstrated, the same outcome follows from a commercially available 
piece of software, COMSOL Multiphysics. An analysis of the origin of the problem, showing that the multival- 
ucdncss of pressure together with its nonintegrability is not inherent in the mathematical formulation and has 
a spurious /numerical nature, is given in Section [5] In Section [6j a method of removing this multivaluedness is 
formulated and numerically verified, first, in the wedge geometry for the steady Stokes flow and then for steady 
and unsteady Navier-Stokes equations in a full-scale dynamic wetting simulation. In Section [7] the results are 
summarized and some their implications are discussed. 



2. Problem formulation 



Consider the two-dimensional steady viscous flow of an incompressible Newtonian fluid, with density p and 
viscosity /i, in a corner confined by straight boundaries located at 9 = and 9 = a of a polar coordinate 
system (r, 6) in the plane of flow and the 'far field' boundary on an arc of a sufficiently large radius r = R. 
Conventionally, we will refer to the 9 = and 9 = a boundaries as the 'solid boundary' and the 'free surface', 
respectively. The 'free surface' will be made genuinely free, with its shape to be determined, in §6.2| 

The flow is driven by the relative motion of a solid at 9 = 0, which slides with speed U parallel to itself, and, 
possibly, also by the far-field conditions. For simplicity, we will assume that the velocity and length scales that 
characterize the flow are such that the Reynolds number Re based on these scales is small. Then as Re — > 0, 
to leading order in Re we may consider the Stokes flow. It should be emphasized that all essential results 
remain valid for the full Navier-Stokes equations since they come from the asymptotic behaviour of the solution 
as r — ► 0. Considering the Stokes flow allows us to demonstrate the method we use to handle the pressure 
multivaluedness more clearly, without additional but nonessential details associated with handling nonlinear 
convective terms. These details are described in §6.2| 

For the problem in question, the non-dimensional Stokes equations for the bulk pressure p and the radial and 
azimuthal components of velocity (it, v) take the form 
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On the solid surface, for a solution not to have multivalued velocity at the corner (see Dussan V & Davis 1974 
Shik hmurzaev|p007 1 , we use the Navier slip condition (see Navier|l823[ ), as opposed to the no-slip condition of 



classical fluid mechanics, and keep intact the impermeability condition for the normal component of velocity to 
the surface: 

du 

d9 



0. 



(0 < r < R, 



0). 



(2.3a,6) 



Here (3 is the dimcnsionlcss 'coefficient of sliding friction' (see Lamb 19321. In the limits — > and (3 — * oo, 
one recovers the conditions of zero tangential stress (free slip) and no-slip, respectively. The value of is 
proportional to a (non-dimensional) 'slip length' that characterizes the region where the velocity field specified 
with the help of the Navier condition (2.3 1) deviates from the velocity field which would have been specified by 
no-slip. 

On the free surface, we have the standard boundary conditions of zero tangential stress and impermeability: 
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0, 



The boundary conditions in the far field can be imposed in different ways. For simplicity, we will make the 
far-field conditions 'passive' and assume that: 

du dv 
dr dr 

This condition is an adaptation for a finite domain of a boundary condition that would specify the asymptotic 
behaviour of t he fl ow field at infinity. As is known (see |Moffatt||1964[ ), this co ndit ion is satisfied if the Navier 
slip condition (2.3 2) is replace by no-slip (i.e. for (3 = 00 in (2.3 2)), so that (2.51 can be seen as a condition 



(r = R, < 9 < a). 



(2.5) 
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that the disturbance caused by finiteness of j3 attenuates in the far field. In computations, for the far-field to 
have a negligible effect on the near-field flow, it is sufficient to put R > 100//3. 
Equations (2.1l-(2.5l fully specify the problem of interest. 



3. Local asymptotics 

The defining feature of our problem is the angle formed by two parts of the boundary and, for future references, 



it is useful to reproduce the leading-order asymptotics for the solution of (2.1 1 (2.5 I as r — *■ (see Shikhmurzacv 



2006 



2007 1 . After introducing the stream function ip by 

r 00 
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equations (2.1 1— (2.2| are reduced to a biharmonic equation A ip — with boundary conditions ( 2.3 1— ( 2.4 1 



taking the form 
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Condition ( |3.2| a) , which is the only inhomogeneous boundary condition in the problem (i.e. the condition that 
drives the flow), suggests looking for the leading-order term of the local asymptotics in the form ip — r 2 F(9), 
which is one of a family of separable solutions to the biharmonic equation of the form ip — r x F(9). After 



substituting ip into the biharmonic equation and boundary conditions (3.2)— (3.3), we have that 

tp = r 2 (B 1 + B 2 9 + B 3 sin 29 + B 4 cos 29) , 
where the constants of integration Bi (i = 1, . . . , 4) are: B\ = — B4 = —(3/4, B2 = —B±/a, B3 = B\ cot 2a. 



(3.4) 



The pressure field obtained from (2.2) using (3.4 1 and (3.1l has the form 



P : 
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(3.5) 



where po is a constant which sets the pressure level. 

The most notable characteristics of the flow are that, as r — *■ 0, the velocity scales linearly with r whilst the 
pressure is logarithmically singular at the corner and does not depend on the angular coordinate 9. These simple 
features provide a quick test to determine whether there are any major issues with a numerical scheme. 



4. Numerical results 



The problem formulated in the previous section was incorporated into a numerical platform based on the 
Galerkin finite clement method. This platform, developed by the authors to tackle a range of microfluidic 
capillary flows, has already been shown to provide novel results in the flow of liquids over chemically patterned 
surfaces (see Sprittles & Shikhmurzaev 20071 12009 1 . As an additional test of the robustness of our numerical 



results (and ubiquity of emerging numerical artifacts) for the flows in the corner regions that we are examining 
here, these results have been verified for test cases using a commercially available code, COMSOL Multiphysics. 
All computations presented below correspond to /? = 10, R= 10; the runs in the process of investigation covered 
a wide range of parameters to ensure that the features described below are invariant with respect to variations 
of these parameters. 

4.1. Details of computations 



The numerical solution of the steady Stokes problem (2.1)-(2.5) uses only a small degree of the platform's 



capability. Given the domain's geometry, it was natural to present the problem formulation in polar coordinates; 
this also made it easy to describe the asymptotic results. However, for a generic free surface flow there is 
nothing to be gained from this implementation and the finite element platform is formulated in the usual 
Cartesian coordinates, with the domain's geometry utilised in the mesh generation (see below). As is relatively 



standard in the finite element approximation of dynamic wetting flows (see Wilson et al. 2006 Lukyanov & 



Shikhmurzaev 20071, we use the P2P1 triangular elements which approximate the velocity quadratically and 
the pressure linearly, thus satisfying the LBB constraint (see|Babuska fc Aziz||1972|. A converged solution to 



the resulting linear equations is obtained after one iteration, using for the inversion the MA41 solver provided 
by the Harwell Subroutine Library. 

The skeleton of the finite element mesh is composed of circular arcs with the corner as their centre, and 
radial-rays which run from r — to r — R. The elements are then tesselated around the skeleton so as to create 
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Figure 1. Streamlines for a = n/4 in increments of ij) = 0.02 with the boundaries at 9 = and at 9 = a corresponding 

to ip = 0. 




Figure 2. Pressure distributions in the vicinity of the corner for a = 7r/4. Left: pressure contours in steps of size 5, as 
the corner is approached with the underlying finite element mesh visible. Right: surface plot of pressure. 



a structured mesh. Moving away from the corner, the distance between arcs is gradually increased by 5% so as 
to provide a high level of resolution near the corner whilst ensuring that the problem remains computationally 
tractable. Due to the nature of the flow it has been found that, as noted in previous investigations of a similar 



nature (see Wilson et al. 2006 Lukyanov et al. 2008 1 , a large number of elements are required in the angular 



direction near the corner in order to accurately capture the dynamics. 

It can be seen in equation (3.51 that the pressure is logarithmically singular as the corner is approached and, 
strictly speaking, one should look to incorporate special 'singular elements' there to capture this behaviour (see 
Wilson et al. 2006 1 . The platform we have developed has the option of incorporating these singular elements 



but, in the results presented here, we have opted against using them in order to simplify our exposition. Their 
implementation has no effect on the conclusions of the forthcoming sections, as verified by our test runs. 



4.2. Acute wedge angles 

First, we consider the case of a < tt/2 and show that our numerical results are in excellent agreement with the 
local asymptotics described earlier. 

The streamlines in Fig.[l]illustrate the general features of the flow: motion is created by the relative movement 
of the solid surface with respect to the corner. The fluid near the solid is pulled out of the corner by the moving 
solid surface and, by continuity, it is replenished there due to the inflow from the far field. 

The pressure distribution near the corner in the form of isobars and, to make it easier to envisage, a 3- 
dimensional plot are shown in Fig. [2] As predicted by equation (3.51 of the local asymptotics, the pressure in 
the vicinity of the corner is independent of 9. 

The quantitative comparison of the computed velocity and pressure with those given by the asymptotics is 
shown in Fig. [3] As one can see, the agreement between numerical and analytic results for the distribution of 
velocity and pressure is excellent. The velocity is linear whilst the pressure, which is plotted in a semi-logarithmic 
frame, is logarithmic with the expected gradient. The pressure constant pg in (3.51 has been chosen to provide 
the best fit, but does not in any way determine the shape or gradient of the curve. 
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Figure 3. Left: comparison of the computed radial velocity u along the liquid-solid interface 1 and liquid-gas interface 
2 with the corresponding asymptotic results 1A and 2A, respectively. Right: comparison of the computed pressure along 
the interfaces to the asymptotic result (dashed line). 




Figure 4. Streamlines for a — 37r/4 in increments of ip = 0.04 with the boundaries at 8 = and at 6 = a 

corresponding to ip = 0. 



4.3. Obtuse angles: multivaluedness of pressure 

In the previous subsection, we have shown that our algorithm gives excellent results for a < ir/2. However, for 
the angles a > ir/2 the situation changes. The same code, as well as COMSOL Multiphysics, that we used for 
comparison, produce results that are markedly mesh-dependent, i.e. the numerical solution cannot be regarded 
as a uniformly valid approximation of the analytic one. 

Fig. [4] shows the picture of streamlines near the corner. The flow is faster than that obtained for acute 
angles but, at first sight, there is no indication of any particular abnormality, at least on a qualitative level. 
However, when we examine the plots of the pressure distribution near the corner shown in Fig. [5] it becomes 
clear immediately that there are severe numerical issues. The smooth ^-independent pressure obtained for acute 
angles is now replaced by two huge spikes of differing signs at the nodes adjacent to the corner. The distributions 
of velocity and pressure along the boundaries near the corner shown in Fig. [6] confirm that we have a serious 
problem. The numerical method implies single- valuedness for both velocity components and the pressure at the 
corner, and computations confirm single- valuedness of velocitjjf] (Fig. [6|. However, as one can see in Fig. |6j the 
velocity differs drastically from what is predicted by the asymptotics, most noticeably by not behaving linearly 
with radius. Furthermore, the asymptotics predicts that near the corner the flow along the free surface is in the 
upstream direction with the radial component of velocity along the free surface being positive. This is in stark 
contrast to the streamline pattern observed in Fig. [4j which shows a regular downstream flow with the velocity 
on the free surface directed towards the corner, as one may intuitively expect. 

Although the problem with the computed velocity distribution is serious enough, the situation with the 
pressure is even more severe. As one can see from Fig. [5] and Fig. [6] the numerical scheme is attempting to 

t Its single- valuedness is ensured by the slip boundary condition (see |Dussan V fc Davis|1974||Dussan|1979 1, whereas, 
numerically, once we assume all functions to be single- valued, the impermeability conditions make the corner a stagnation 
point. 
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Figure 6. Left: comparison of the computed radial velocity it along the liquid-solid interface 1 and liquid-gas interface 
2 with the corresponding asymptotic results 1A and 2A, respectively. Right: computed pressure along the interfaces. 



approximate a solution which is both singular and multivalued at the corner. As one approaches the corner 
along different radii, the pressure behaves differently and tends to plus infinity along the radii closer to the 
free surface and to minus infinity along those closer to the solid boundary. Since the numerical scheme imposes 
(artificial) single-valuedness at the corner by stating that pressure has one value at the corner node, the code 
tries to reconcile the calculated pressure with this requirement over the first row of elements, thus creating 
the cliffs shown in Fig. [5] This makes the numerical results mesh-dependent: the smaller one makes the size of 
the first row of elements, the nearer the multivalued solution approaches the corner, and the higher the cliffs 
of pressure become. Obviously, the mesh-dependent numerical results cannot be regarded as a uniformly valid 
approximation of the actual solution to (2.1 1 — ( 2.5 1 and hence have to be rejected. 

Furthermore, as shown by the semi-logarithmic plot in Fig. [6] the singularity of pressure at the corner 
is stronger than logarithmic (by plotting the pressure in the log-log frame one can show that it is actually 
algebraic with the exponent dependent on a). Besides being at odds with the asymptotics, such behaviour 
of pressure poses a fundamental problem for computations of free-surface flows with moving contact lines for 
which the corner flow considered here is but an element: a nonintegrably singular pressure distribution along 
the free surface makes corrections to the wedge shape divergent and hence does not allow one to compute the 
free surface shape at all (see |Shikhmurzaev||2007 1 . 

It should be emphasized that the difficulties we are describing are above the level of a particular numerical 
implementation of a particular algorithm. It is not only that the code we used has been thoroughly validated 
(see the previous section); additionally, the commercially available code COMSOL Multiphysics has also been 
applied to the wedge flow problem for both acute and obtuse angles, and in all cases identical results have been 
recorded. 

The encountered problem turns out to be resilient to standard approaches which have been described in the 
literature as successful remedies to spurious numerical effects. The first approach to be tried is to incorporate 
the asymptotic results described earlier in the numerical scheme to 'come out' of the corner, thus bypassing 
any spurious effects that might have been caused by the geometry. The asymptotics can be matched with the 
numerical solution outside a given radius in a variety of ways. This approach has proved successful in similar 



situations in which corner singularities exist (see Shi et at. 2004 1. However, it has been found that for 
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problem such an alteration of the scheme merely shifts the pressure cliff to the arc at the radius where the 
asymptotics has been applied. This is of no help whatsoever, since the results have all the deficiencies listed 
above. 

Another standard approach is to impose penalties of various forms, similar to those used with equal-order 
interpolation to circumvent the LBB condition (see Baer et «L|2000 ), but such a method simply fails to drive the 
code into a mesh-independent solution. The extent to which the penalties flatten the pressure cliffs is exactly the 
extent to which the numerical solution (to the 'penalized problem') departs from that which would satisfy the 
original problem discretized using the standard FEM. In other words, instead of driving the code to the 'right' 
solution, the penalty becomes part of the differential equations the code is solving, i.e. the method effectively 
replaces the original problem with a different one. 

The spurious pressure behaviour which we have shown throughout the paper has previously been observed 
and designated the subject of future research in a paper on the finite element simulation of curtain coating 
(see Wilson et a?.||200l" l. As is pointed out in this paper, a reason that the problem may have not been treated 
in other investigations is a lack of computational resolution. A rough indication of the region which must be 
resolved by any numerical scheme is the slip length, given by l/f3. Tests show that at least 100 nodes should be 
used along each radial-ray within this region. With the graded mesh we have used and for modern computers 
this requirement causes no problem. For a mesh with a single element size, achieving this is more difficult, as 
the numerical cost of having every element of size 0.01//3 could render the computational task untractable. It 
may be that this is the case in a number of older publications in which computational power proved a major 
problem (see Kistler & Scriven 1983 Christodoulou & Scriven 19891. In our calculations presented here, the 
smallest element has the size of 4 x 10 _y , and we have 206 nodes inside the slip length along each radial-ray. 



5. Origin of the pressure multivaluedness 

At this stage, the question is whether the difficulty, albeit resilient to conventional alterations of the algorithm, 
is a numerical artifact, or whether at a > ir/2 the computed features are a genuine property of the mathematical 
problem. The robustness of the velocity distribution and the fact that, contrary to the asymptotics, its radial 
dependence is strongly nonlinear, with a singular radial derivative at the corner, suggest that here we might 
be dealing with an eigensolution, i.e. a solution satisfying homogeneous boundary conditions, v = du/d9 = 
(9 = 0,a), that superimposes on the one described by the asymptotics. Then, if the eigensolution for the 
velocity has indeed singular radial derivatives at the corner, the computed pressure behaviour could result from 
numerical errors in calculating the velocity field corresponding to this eigensolution. 



An eigensolution to our problem has been known for a long time (see Moffatt 1964 P. 14). In terms of the 



stream function, it has the same separable form as the asymptotic solution described in Section [3] 

V^ e = Ar x sin(X9), A = n/a, (5.1) 

where A is an arbitrary constant. 

This solution is promising for two reasons. Firstly, it is only for a > ir/2 that we have A < 2 and hence this 



solution can dominate the one described by the local asymptotics (3.4) in the near field. Secondly, it exhibits 
the numerically computed non-linear radial dependence of velocity. Note that, although the velocity of the 
eigensolution tends to zero at the corner, its radial dependence scales like r A_1 , and hence for it/2 < a < ir the 
derivatives of velocity in the radial direction are singular at the corner. 



Crucially, the eigensolution (5.1l produces a globally constant pressure. This simplicity allows an ideal op- 
portunity to check whether the nonintegrability and multivaluedness of pressure computed earlier are indeed 
numerical artifacts. In order to do this, we consider the flow in a corner region with the Navier slip condition 
replaced with zero tangential stress, i.e. with du/d9 — and v — on both boundaries 9 = and 9 = a. The 
flow can then be generated only by the boundary conditions in the far field which we set using the eigensolution 



(5.1 1, where, for simplicity, we use A = 1/A: 

u = r A_1 cos(A0), v = -r A ~ 1 sin(A0) (r = R, < 9 < a). (5.2) 

Then the eigensolution is the exact global solution to our test problem. As for the pressure distribution, once 
the pressure level has been set, say, to zero, then one will have p = in the whole domain. 

Once this problem is computed numerically, the situation becomes clear. Fig. [7] shows the isobars and a 
3-dimensional plot obtained using COMSOL Multiphysics, again with the P2P1 element. The isobars and the 
3-dimensional plot resemble those which we have already seen for the Navier condition in the previous section, 
with the same cliffs in the pressure profile at the nodes nearest to the corner. This is a remarkable result given 
that we know the global solution is in fact p = 0! The results obtained using our numerical platform and those 
of COMSOL are in perfect quantitative agreement in the far field and are similar in the near field; there is little 
point in comparing specific values in the near field as the results are seen to be mesh dependent for both codes. 
An interested reader can easily reproduce the results using COMSOL, where the condition of impermeability 
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Figure 7. Pressure distributions in the vicinity of the corner for a = 37r/4 as calculated by COMSOL Multiphysics. 
The underlying finite element mesh is visible in both plots. Left: pressure contours as the corner is approached. Right: 
surface plot of pressure. 



and zero tangential stress is selected as a boundary condition by choosing boundary 'wall' and then choosing 
the type 'slip'. 

The underlying mesh, visible for both plots in Fig. [7J has been generated by COMSOL using adaptive grid 
refinement techniques and it is unstructured. The same spurious features of the numerical solution have been 
encountered for both the unstructured mesh of COMSOL and the structured mesh of our numerical code, thus 
suggesting that any alternative mesh design would not allow a route out of the problem. Additionally, COMSOL 
offers the choice of ten different elements, and, as we have tested, the choice of any of these elements does not 
affect the overall picture. 

Therefore, we may now conclude that the pressure multivaluedness is spurious and has a numerical origin. In 
the next section, we will show how the information we have about the cigcnsolution as the underlying cause of 
the problem can be used to resolve it. 



6. Removal of pressure multivaluedness 

First, we will describe the method of removing the pressure multivaluedness from the numerical solution of 
the problem formulated in Section 2, i.e. for the Stokes flow in a wedge region, and then show how the method 
can be 'localized' to apply to general 2-dimensional Navier-Stokes flows, both steady and unsteady, where an 
angle formed with the solid surface by a priori unknown free boundaries is but one element. 

The key idea of the method is very simple. In the situation where, as in our case, the eigensolution is the 
cause of the pressure multivaluedness, we can subtract this solution from the problem and use the degree of 



freedom it offers, i.e. arbitrariness of A in (5.1l, to ensure single- valuedness of pressure at the corner. Once the 



supplementary to the eigensolution velocity and the pressure fields are computed, we can put the cigcnsolution, 



i.e. the velocity field described by (5.1 1 and uniformly zero pressure, back in. The resulting solution will have the 



analytically known eigencomponcnt of the velocity field superimposed on the computed 'supplementary' flow 
and the computed single-valued pressure. The resultant combined solution will satisfy the original equations 
and boundary conditions. 

6.1. Simplest case: Stokes flow in a corner region 
For the problem formulated in Section [2j consider the velocity and pressure as sums of the eigensolution 

u e = AAr A_1 cos(A0), v e = -AAr A_1 sin(A0), p e = 0, (A = n/a), (6Aa,b,c) 

and the components to be computed (hereafter these are marked with a tilde): 

(u,v,p) = (u ei v e ,p e ) + (u,v,p). (6.2) 



The constant A in ( |6.1| a,6) is yet to be specified. 

Since the eigensolution satisfies the Stokes equations (2.1)— (2.2| and the free-surface boundary conditions 



(2.4 1 exactly, one has that u, v and p have to satisfy the unaltered Stokes equations (2.1 1— ( 2.2 ) and boundary 



conditions (|2.4|), whereas the boundary conditions on the solid surface and in the far field for these variables 



will take the form: 



09 



= r/3 (AXr*- 1 + u-l), v = (0 < r < R, 9 = 0), 



(6.3a,6) 



AX(X — l)r sin(A6>), (r = R, < 6 < a). (6.4a, 6) 



9 £ = -A\(X-iy->co S (A6), g 
To complete the problem formulation for u, v and p, we must add an equation to account for the additional 
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Figure 8. Top: the computed streamlines in increments of ip = —0.04. Bottom: the streamlines with the eigensolution 
superimposed to produce the full solution with increments of ip = 0.04. In both plots a — 3n/4 with the boundaries at 
9 — and at 8 = a corresponding to ip = 0. The constant A — 1.3026. 



unknown constant A. To do so, we impose a condition that the pressure is single valued at the corner: 

lim H = 0. 



(6.5) 



Qualitatively, this condition can be explained as follows. The method is based on taking the eigensolution out 
of the total and hence ensuring that (m, v) do not have the singularity of the radial derivatives at the corner, 
and the numerical error in their computation will not give rise to the errors in computations of the pressure 



resulting in its multivaluedness. By imposing (6.5 1, we are effectively ensuring that the eigensolution is taken out 



fully from the viewpoint of what this subtraction is aimed at achieving. In other words, out of a one-parametric 



family of cigcnsolutions, parameterised by the constant A, condition (6.51 selects the one that underpins the 
flow we are considering. 



A simple way for us to impose (6.5 1 numerically is to demand that the pressures at the nearest to the corner 



nodes on the free surface and on the solid boundary are equal, as in our mesh they are equidistant from the 
corner. 



Equations ( 2. 1 )-( 2.2 ) , ( 2.4 ) , (6.31(6.5) have been solved using our numerical platform with the same solution 



procedure as before. The streamlines of the supplementary flow obtained from the computation are shown in 
Fig. [8] alongside the streamlines to our original problem formulated in Section 2, obtained using (6.2 1. Although 
the underlying asymptotic solution predicts that there must be flow reversal near the corner, which is replicated 
in our numerical solution for (u,v), this feature is blown away by the strength of the eigensolution when it is 
superimposed on top. 

The pressure plots in Fig. [9] show a qualitative transformation from those observed in previous sections: the 
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Figure 9. Pressure distributions in the vicinity of the corner for a — 37r/4 using our numerical platform with the 
eigensolution removed and A — 1.3026. Left: pressure contours in steps of size 2.5 as the corner is approached, with the 
underlying finite element mesh visible. Right: surface plot of pressure. 




Figure 10. Left: comparison of the computed radial velocity u along the liquid-solid interface 1 and liquid-gas interface 
2 with the asymptotic results 1A and 2A, respectively. Right: comparison of the computed/actual pressure along the 
interfaces to the asymptotic result (dashed line). 



pressure is now single valued, smooth and exhibits no mesh-dependence. Both the isobars and the 3D 'surface 
plot' are now similar to those obtained earlier for acute angles (Fig. [2]). 

The comparison with the asymptotics of Section 3 is shown in Fig. |10| The agreement of pressure is visibly 
excellent. The fit between the asymptotics and the radial velocity along the interfaces, u, is now confined to a 
smaller region than that previously observed for acute angles. This is no surprise, given that the eigensolution 



now enters the Navier conditions (6.3 2), where the first term in brackets on the right-hand side becomes small 
compared to the last term only very close to the corner. 

Interestingly, the dominance of the eigensolution in the combined solution for obtuse angles ensures that the 
velocity field near the corner is almost anti-symmetric about the centre line 9 — a/2, whereas for acute angles 
this is not the case. 

Thus, for the present problem of flow in a corner, we have fully resolved the situation. However, when consider- 
ing more complicated flows which involve a corner only as an element , the method of removing the eigensolution 
from the problem formulation throughout the entire domain could become unnecessarily complicated; the eigen- 
solution is only important near one corner and yet it will make artificial contributions to equations and boundary 
conditions in the whole domain. A more reasonable approach would be to design a 'local' variant of the method 
to remove the eigensolution only near the corner, where its presence creates the unwanted numerical artifacts, 
leaving the rest of the flow domain intact. This variant of the method is considered below. 

6.2. General case: Navier-Stokes equations in domains with curvilinear free boundaries 

We will begin by describing the pressure regularization method in its generic form and then illustrate it by 
considering two test problems: (a) a steady propagation of a meniscus in a channel with plane-parallel walls 
and (b) an unsteady spreading of a (2-dimensional) liquid drop over a solid surface. In both cases, we will be 
describing only the details related to the implementation of the method and leave out the standard elements of 
the problem formulation for these free-surface flows. 

Now, we have that the bulk flow is described by the full Navier-Stokes equations up to the corner (which, in 
accordance with the literature, we will now call the 'contact line') and the position of the curved free surface is 
a priori unknown. 

In order to 'localize' the pressure regularization method, we now introduce an artificial 'internal boundary' 
that will separate the 'inner' region near the contact line where the regularization procedure will be used from 
the rest of the flow domain (the 'outer region'). Inside the inner region, we remove the eigensolution prior to 
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computing the supplementary components of the velocity and pressure, and then superimpose it back, as we did 
in the previous subsection: this prevents the spurious pressure behaviour described in Sections 4 and 5. In the 
outer region, there is no need to alter the standard numerical approach to solving the appropriate equations. 
Solutions in both the inner and the outer region are computed simultaneously and have to be matched at the 
internal boundary. 

Thus, after decomposing the solution in the inner region into the sum of the eigensolution and the supple- 
mentary part to be computed, 

(u,p) = Oe,0) + (u,p), (6.6) 

we need to: 



(i) Take into account, as we did in [6.1 the contribution of the eigensolution to the boundary conditions for 
u on the solid surface. 

(ii) Take into account the contribution of the eigensolution to the bulk equations for u and p in the inner 
region arising from the fact that the eigensolution satisfies the Stokes but not the Navier-Stokes equations. 

(iii) Take into account the contribution of the eigensolution to the boundary conditions for u and p on the 
free surface that appears due to the fact that now the free surface, now being genuinely free, is not necessarily 
planar. 

(iv) Formulate the matching conditions at the internal boundary that would link the ii and p with the outer 
flow. These conditions are necessary to calculate solutions in both the inner region and the outer region; these 
calculations are carried out simultaneously 

After non-dimcnsionalizing the problem using characteristic length L and velocity U scales from the global 
flow, on the flat solid surface, which moves at non-dimensional speed U w , we have the same Navier-slip and 
impermeability conditions as in |6.1| 



n ■ P ■ (I nn) = (3(u- 



U w ) ■ (I — nn) , (u + u e — U w ) ■ n = 0, 



(6.7a, b) 



where P = —pl+ Vit + (Vit) T is the stress tensor of the supplementary solution; J is the metric tensor, n is 
the internal normal to the interface, so that the tensor (J — nn) extracts the tangential component of a vector. 

(6.8) 



In (6.7 1), we have used that for the stress tensor of the eigensolution, 



Pe 



V« e + (V«e) 



one has n ■ P e ■ (I — nn) = 0. 

The bulk equations for u and p to be solved in the inner region take the form 



V-tt = 0, 



Re 



d(ii + u e ) 
dt 



+ (u + u e ) ■ V (it, + u e ) 



= V -P + Stg, 



(6.9a, 6) 



where St = pL 2 g/(^U) is the Stokes number based on the acceleration due to gravity g and g is a unit vector 
in the direction of the gravitational force. In equations ( |6.9[ 3,6), we have used that the eigensolution satisfies 
the Stokes equations, i.e. V • u e = V • P e = 0. 

On the free surface, which is defined by f(x,t) = 0, where x is the position vector and the function / is to 
be found, we have the standard dynamic and kinematic conditions: 



Can - (P + P e 



nV 



dl 

dt 



+ (w + O-V/ = 0, 



(6.10a, b) 



where Ca — fiU/cr is the capillary number based on the constant surface tension coefficient a of the liquid-gas 
interface. 

In the numerical implementation, care must be taken in the evaluation of the term n ■ P e in (6.10 2) as, 
although it is intcgrablc, it is singular in the limit r — »• 0. The term is best evaluated by calculating the stress 
tensor analytically, rather than using the finite clement approximation, or whatever other discretization has 



been chosen, for the derivatives of the cigcnsolution's velocity components in (6.8 1 



An internal boundary separating the inner region from the outer flow should lie sufficiently far away from the 
contact line for the eigensolution to be well within the inner region and at the same time not too far for the 
regularization method to be localized, as opposed to applied to an unnecessarily large region of the overall flow. 
Another consideration is how the position of the internal boundary correlates with the computational mesh. In 
our numerical method this is most easily achieved by using one of the arcs formed by the edges of the finite 
elements; the method is equally applicable for an algorithm with an unstructured mesh, though the ease of 
defining the internal boundary will be lost. 

At the internal boundary, we enforce continuity on the velocity and the stress: 



u f 



u = u n 



n; 



P) = m-P a 



(6.11a,6) 
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Figure 11. Streamlines near the contact line of a propagating meniscus. The free surface and solid surface both cor- 
respond to ip = 0. A is computed to be 1.4551. Top: streamlines of the calculated velocity field with increments of 
ip — —0.02 in the inner region and if) = 0.02 in the outer region. Bottom: plot showing the actual streamlines obtained 
after superimposing the eigensolution with the calculated solution in the inner region. 



where n.^ is a normal to the internal boundary and the subscript out marks the velocity and stress in the outer 
region that also have to be computed. 

The ability of the finite element method to incorporate stress boundary conditions in a natural manner 
ensures that the conditions of continuity on the actual solution (6.111 can be satisfied without dropping any 
other equations. 

In order to illustrate how the method works, we consider two model two-dimensional simulations using a 
Cartesian frame x — (x, y) with a solid surface at y = 0. In our simulations, we fix the parameters to Ca = 0.1, 
Re = 1, St = 1 and /3 = 10. As our interest here lies in the numerical approximation of these flows, as opposed 
to a detailed comparison with experiment, we may consider, for simplicity, the dynamic contact angle, which 
was previously labelled a, to have a fixed value of 9d = 37r/4. 

The first problem is the steady motion of a meniscus that a fluid-gas interface forms between plane-parallel 
walls. In a frame moving with the contact line the problem is time-independent and hence the time derivatives 



in (6.9 b) and (6.10 b) axe zero. The velocity of the solid substrate in the moving reference frame is U w = (1,0). 
Fig. 11 shows the velocity fields as they have been computed and the resulting (combined) velocity field. In 



the top picture, we can see the outer flow and a much weaker supplementary flow in the inner region. When the 
eigensolution is superimposed back on top of the supplementary velocity field in the inner region, the streamlines 
are seen to not feel the presence of the internal boundary — the matching conditions work perfectly leaving no 
'scar' on the flow. As one can observe, peculiarities of the underlying (supplementary) flow in the inner region, 
such as flow reversal, are of little consequence as their effect is negligible compared to that of the eigensolution. 
Finally, we show that our method is equally applicable to time-dependent flows. As an illustration, we consider 



Numerics of Viscous Flows in Domains with Corners 



13 




0.0 0.5 1.0 

Figure 12. Snapshot at t — 1.5 of the streamlines in a spreading liquid drop. The solid surface corresponds to ip — 0. A is 
computed to be 1.6025. Streamlines in all plots are in increments of ip — —0.05. Top: plot showing the computed velocity 
field in both the inner and outer regions. Bottom: plot showing the actual streamlines obtained after superimposing the 
eigensolution with the calculated solution in the inner region. 

the spreading of a two-dimensional liquid drop over the solid surface, with the axis of symmetry at x = 0. The 
drop is driven from its initially cylindrical shape by gravity. In Fig. |12| we show a snapshot of the streamlines 
near the contact line in a frame moving with the contact line. The free surface is in a state of evolution and hence 
no longer represents a streamline. Again, a comparison of the two plots in the figure show that the position of 
the transition line does not effect the overall flow. 



7. Conclusion 

We have shown that straightforward application of a standard numerical method to a seemingly ordinary 
problem can lead to completely unacceptable numerical artifacts, despite the fact that the conventional prelim- 
inary asymptotic analysis of a possible source of difficulties (the corner formed by smooth parts of the domain's 
boundary) does not flag up any concerns. 

A surprising result from the present study is that errors in approximating the velocity field manifest themselves 
as spurious behaviour in the pressure field. It is the presence of an underlying eigensolution which creates these 
numerical artifacts in the pressure distribution, despite the fact the eigensolution itself has a globally constant 
pressure. Therefore, preliminary analytic work should always include a search for possible cigcnsolutions, an 
analysis of their properties and an attempt to numerically approximate them individually. Only using this 
procedure have we been able to determine the origin, and hence solution to, the problems presented in Section]?] 

We have shown that, if the eigensolution is not removed prior to computations, one invariably ends up 
with huge pressure spikes whose position and magnitude are both mesh-dependent. The numerical analysis 
suggests that the cause of this numerical instability is the errors in approximating the velocity gradient, as 
in the eigensolution this gradient is singular at the corner. The 'mechanism' of this error generation and the 
subsequent instability require a thorough mathematical investigation, which lies outside the goal of the present 
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paper. Such an investigation would, hopefully, point out other situations where standard numerical methods 
could run into trouble. 

The developed method of removing numerical artifacts in the pressure distribution is not only successful 
with respect to the model case of a steady two-dimensional Stokes flow in a corner region, but, as is shown, it 
admits a straightforward generalization which makes it applicable to a general case of unsteady free-boundary 
Navier-Stokes flow. In practical applications to problems of dynamic wetting one often has a situation where 
in the process of computations the contact angle varies in a wide range, from the angles where a standard 
numerical code produces no artifacts to those where the pressure spikes and multivaluedness invariably appear. 
General-purpose numerical algorithms should be developed so that the present method is turned off for acute 
angles, where the pressure is naturally single-valued, and switched on for obtuse angles to suppress spurious 
numerical behaviour. 

The failure of standard numerical algorithms to approximate flow in a corner is not limited to the formulation 
considered in the present paper. For example, an alternative mathematical approach to the moving contact-line 
problem (see Dussan V 1976 Zhou & Sheng 1990 Somalinga & Bosc 2000) is to prescribe the velocity as an 
explicit function of distance from the corner, so that, in a frame moving with the contact line, the velocity at 
the contact line is zero and it tends to the speed of the solid in the far field. When the Dirichlet conditions of 
this type are applied instead of the Navier-slip (i.e. Robin-type) condition that we have examined, one observes 
spurious numerical features, though in a different range of the corner angles to those in the present paper. This 
case is examined in a forthcoming paper. 
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